Negative Temperature States in the Discrete Nonlinear Schrodinger Equation 
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We explore the statistical behavior of the discrete nonlinear Schrodinger equation. We find a 
parameter region where the system evolves towards a state characterized by a finite density of 
breathers and a negative temperature. Such a state is metastable but the convergence to equilibrium 
occurs on astronomical time scales and becomes increasingly slower as a result of a coarsening 
processes. Stationary negative-temperature states can be experimentally generated via boundary 
dissipation or from free expansions of wave packets initially at positive temperature equilibrium. 
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The discrete nonlinear Schrodinger (DNLS) equation 
has been widely investigated as a scmiclassical model 
for Bosc-Einstcin condensates (BEC) in optical lattices 
and for light propagation in arrays of optical waveguides 
[ij. It exhibits peculiar solutions in the form of lattice 
solitons, also known as breathers, or intrinsically local- 
ized nonlinear excitations [2[. In the first statistical- 
mechanics study of the DNLS equation, Rasmussen et 
al. [3J] identified a region R n in the parameter space that 
was conjectured to be characterized by negative temper- 
ature (NT) states. NT states have attracted the curiosity 
of researchers sincepioneering works in systems of quan- 
tum nuclear spins [4|]. In a series of recent papers, Rumpf 
argued that there are no NT equilibrium states in R n . In 
fact, he showed that the maximum entropy (i.e. equilib- 
rium) state is characterized by an infinite-temperature 
background superposed to a single breather that collects 
the "excess" energy [5[ . Convergence to such equilibrium 
states, however, has never been observed, since the tran- 
sient may last for astronomical times [lfj . Therefore it is 
legitimate to ask what happens over physically accessible 
time scales. The occurrence of long transients is not rare 
in statistical mechanics: it may be due to coarsening, nu- 
cleation, presence of free-energy barriers, or the stability 
of some modes. Within the physical setups that are clos- 
est to the DNLS, a nonexponcntial relaxation was first 
found in chains of nonlinear oscillators and shown to orig- 
inate from the presence of long-lived localized solutions 
[6|| . Slow relaxations have been also found in Heisen- 
berg spin chains and traced back to the existence of two 
conservation laws [7J . This latter example is particularly 
instructive, since it is the context where theoretical argu- 
ments, later applied to the DNLS, have been developed. 
The approach is well suited for the identification of the 
equilibrium state by maximizing the entropy given the 
constraints imposed by the conservation laws. It is how- 
ever of little help to investigate the convergence proper- 
ties to equilibrium. In this respect, the implementation 
of a microcanonical Monte Carlo (MMC) method turns 



out to be quite useful, in that it allows uncovering (in 
the so-called weak-coupling limit) a first source of "slow- 
ness" in an underlying coarsening phenomenon that fol- 
lows from the subdiffusive behavior of the breather am- 
plitude. As a result, we find that the density of breathers 
vanishes in time as t~ a with an exponent a » 0.37 that 
is reminiscent of the exponent 1/3 of the Cahn-Hilliard 
model [8|. Note, however, that dynamical mechanisms 
which strongly inhibit the diffusion of the breather am- 
plitude are even more effective. With the help of an ex- 
plicit expression for the microcanonical temperature [9( , 
we are able to show that a broad class of initial con- 
ditions (IC) of the DNLS equation converges towards 
a well defined thermodynamic state characterized by a 
negative temperature and a finite density of breathers. 
There is no reason to dismiss the theoretical arguments 
of [5| , although it appears that the dynamical freezing of 
the high-amplitude breathers slows down the evolution 
so much to make the convergence to equilibrium unob- 
servablc. Altogether this phenomenon is reminiscent of 
aging in glasses, but a more detailed analysis is required 
to frame the analogy on more firm grounds. Finally, we 
briefly comment on a simple strategy to generate exper- 
imentally metastable NT states in BEC in optical lat- 
tices and in arrays of optical waveguides. Our scheme 
is far simpler than the thermalization methods recently 
proposed for BEC trapped in lattices which require non- 
trivial Mott states followed by sweeps across Feshbach 
resonances to invert the sign of the atomic interaction 

EH- 

Model and state of the art. With reference to the stan- 
dard canonical coordinates, the Hamiltonian of a DNLS 
chain writes 

H(p, ?) = ]T(p 2 + q)f + 2 J^(pjP j+1 + q 3 q 3+1 ) , (1) 
J J 

where the index j numbers the sites from 1 to TV and pe- 
riodic boundary conditions are assumed. The sign of the 
first term in the r.h.s. is assumed to be positive, since 
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FIG. 1: (Color online) Phase diagram (a,h) [jg. We show the 
isothermal lines T = (solid blue), T — 1 (solid purple) and 
T = oo (dashed red). In the inset: average trajectory (over 
an ensemble of 500 realizations) of a chain with N = 200 in 
the presence of boundary dissipation. 



FIG. 2: (Color online). Average breather density p versus 
time t in a DNLS chain as from MMC simulations in the 
weak-coupling limit. Circles, squares and triangles refer to 
TV = 400, 1600, and 6400, respectively. The straight line with 
slope 0.37 is the result of a power-law fit. 



we refer to the case of a repulsive-atom BEC, while the 
sign of the hopping term is irrelevant, due to the symme- 
try associated with the canonical (gauge) transformation 
</>j — > (f>j +7T (where pj + iqj = UjC 1 ^ represents the wave 
function at site j). The model possesses two conserved 
quantities: the energy H and the number of particles (or, 
wave action) A = Y]- \A- The thermodynamic properties 
of this model are summarized in a phase-diagram of en- 
ergy (h) versus particle (a) density per site [3j. In Fig. [I] 
one can identify three regions: (i) a forbidden region Rf, 
that lies below the ground state; (ii) an intermediate pos- 
itive temperature R p region; (iii) the above mentioned 
region R n , that is investigated in this Letter. The pres- 
ence of low relaxation phenomena in R n was first noticed 
in [3j by monitoring the probability density of the local 
amplitudes Aj = u?. The reason was traced back to the 
presence of breathers (that arise as a result of a modula- 
tional instability J7J, [l3[) and attributed to the weakness 
of their coupling to the background. 

Microcanonical Monte Carlo simulations. We start our 
analysis of DNLS dynamics from the weak-coupling limit 
(i.e. large particle densities). In this regime, the region 
R n corresponds to h/a 2 > 2 [5| (with our normaliza- 
tions). The motivation for a MMC study is a direct inves- 
tigation of entropic effects during the time evolution. In 
order to mimic the deterministic evolution, we have intro- 
duced a local updating rule which applies to a randomly 
selected triplet of consecutive sites. Since it is necessary 
to leave the number of particles Aj = A7-1 + Aj + Aj+\ 



and the energy H? = A~ 



4 



A J+1 unchanged, the 



new configuration must lie along the intersection between 
a plane and the surface of a sphere, while restricted to 
the positive octant (Aj > 0). Depending on the rela- 
tive amplitude of H? and Ap this line may be either a 
full circle or made of three separate arcs. In the former 
case, we select a point on the circle (assuming uniform 
probability); in the latter one, we select a point within 



the same arc as of the IC. It is easy to verify that de- 
tailed balance is satisfied. After having observed that 
the above algorithm yields a fast convergence towards 
the expected equilibrium state in the T > region, we 
have performed some simulations in R n . The IC is fixed 
by assigning a random amplitude Aj (with a flat dis- 
tribution between and 1) to all sites and adding 40/3 
units to 5% of sites randomly chosen over the chain. This 
choice implies that h/a 2 ps 7 > 2, i.e. the IC lies inside 
R n . In Fig. [2] we plot the MMC evolution of the frac- 
tion p of sites, where the amplitude Aj is larger than the 
threshold r) = 10 [14[ for three different system sizes [15| . 
The three curves nicely overlap and show for large time 
t a clear power-law decay p ~ t~ a , with an exponent 
a rj 0.37. Time t is given by the number of MMC steps 
divided by N. This means that the system converges, 
although slowly, to the equilibrium state characterized 
by a single breather (i.e. p = in the thermodynamic 
limit). In order to gain further insight, we have looked 
at the MMC evolution in two different ways. In Fig. (3Ji, 
we plot the position of the breathers of amplitude larger 
than 77. The breathers are basically static and either 
terminate abruptly, or eventually exhibit a weak mutual 
attraction. The representation in panel Fig. [3Jd is more 
illuminating. There, we plot the sum of the "excess" 
energies H* = YH=i(Af - r] 2 )0(A 2 - if) (8 is the Heav- 
iside function) in correspondence of each breather. The 
abscissa of the rightmost breather is the total excess en- 
ergy contained in the breathers. After an initial transient 
during which the excess energy is constant, the breathers 
perform a sort of Brownian motion in energy space un- 
til they collide and merge. This (sub)diffusive motion 
reveals a slow convergence that is due to a a likewise 
coarsening phenomenon. It is intriguing to see that the 
value of a (0.37) is relatively close to the critical expo- 
nent of the Cahn-Hilliard equation (1/3). We leave the 
task of identifying the correct universality class to future 




FIG. 3: (a) Space time representation of MMC dynamics of 
breathers in the NT region, (b) Same as in panel a, with the 
label j replaced by Hj, see text. 
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FIG. 4: Evolution of the local amplitude in an NT states (dots 
correspond to points where Aj > 10) for a = 1 and h = 2.4. 



investigations. 

DNLS Dynamics. The MMC study has shown that 
in the weak-coupling limit, the convergence is character- 
ized by a coarsening process: the fewer the breathers, the 
weaker their interaction and the slower the convergence 
process. The uncoupled limit is, however, rather pecu- 
liar, as it lacks the hopping mechanism. In the natural 
case of a finite coupling, we have turned our attention 
to simulations of the deterministic DNLS since effective 
implementations of MMC rules are not feasible. Fig. |4] 
presents the space-time evolution for a chain of 4096 sites 
for a = 1 and h = 2.4. To emphasize the breather dy- 
namics, we plot only the lattice points, where the instan- 
taneous mass Aj is larger than i] = 10. Thanks to a sym- 
plectic 4th-order algorithm of Yoshida type [16[ , we have 
simulated the dynamics over time scales much longer 
than in the previous numerical studies. In Fig. @] one 
can see that, like in the MMC simulations, the breathers 
do not basically move; at variance with the uncoupled 
case, however, we observe spontaneous birth and death 
of breathers as in a standard stationary process. This 
is due to the presence of a finite interaction (hopping) 
energy: the background can store excess energy in the 
phase differences of neighbouring sites and this implies 
that breathers can spontaneously nucleate. In order to 
better investigate the convergence properties of the dy- 
namics, we have built four different sets of ICs, all with 
the same energy and number of particles but substan- 
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FIG. 5: (Color online) Evolution for a = 1, h = 2.4 of (a) the 
breather density p and (b) the inverse temperature /3 = 1/T 
starting from different ICs for a chain with N = 4096: solid, 
long-dashed and dashed lines refer to a suitable homogeneous 
background plus 0, 4 (amplitude 20) and 8 (amplitude 16) 
breathers, respectively, while the dot-dashed curves refer to 
half a chain empty. Both p and /3 are averaged over a moving 
window of 10 time units. Notice that in panel (b) the two 
dot-dashed lines correspond to spatial averages in the initially 
empty and filled sectors of the chain, (c) Breather probability 
density after a transient of 5 • 10 6 time units. Data are not 
collected below the threshold r\ = 10. In the three panels 
the dotted lines refer to homogeneous ICs and N — 2048. 
All data sets have been also averaged over a dozen different 
realizations of each IC. 



tially different macroscopic structures (see the caption 
of Fig. \5§. First we have monitored the evolution of 
the density of breathers (identified by setting 77 = 10). 
The resulting average results for TV = 4096 are plotted 
in Fig. [5^,. At variance with the simulations described in 
Fig-Hi the density appears to convergence towards a com- 
mon finite value, i.e. towards a multi-breather stationary 
state. Moreover, the dotted line corresponding to a simu- 
lation with N = 2048 indicates that the asymptotic value 
of p is independent of the system size (and approximately 
equal to one breather per thousand sites). We have also 
monitored the temperature. In this model, the standard 
kinetic definition does not apply, since the Hamiltonian 
is not decomposable into kinetic and potential energy. 
Nevertheless, one can use the microcanonical definition, 
T _1 = dS/dH, where S is the thermodynamic entropy 
and the partial derivative must be computed taking into 
account the existence of two conserved quantities. This 
task requires long and careful calculations with a final 
non-local expression [9J. We have preliminarily verified 
that the definition of [9j works for positive temperatures 
when applied to the full and to short sub chains. The 
results in R n are plotted in Fig. [5t>, where one can see 
that j3 = 1/T converges to a common value for all of the 
four different classes of ICs (either from above or from 
below), with negligible finite-size effects (once again the 



dotted curve refers to a chain of half length, N = 2048). 
Altogether, one can conclude that on the time scales that 
are numerically accessible (of order 10 7 ), the DNLS con- 
verges towards a multi-breather state characterized by 
a finite density of breathers and a well-defined negative 
temperature. 

In practice, what happens is that the breathers are 
prevented from becoming "too large". In Fig. [5b, one 
can indeed see that the pseudo-stationary distribution of 
breather amplitudes is effectively localized below a = 22 
and, more importantly, the distribution is basically inde- 
pendent of the system size. Altogether, this means that 
there is a dynamical process which "screens" the high 
amplitude region. One can argue that this dynamical 
effect is due to the low efficiency of the energy-transfer 
among breathers interacting through the background, a 
manifestation of their intrinsic dynamical stability. As a 
consequence, the evolution is "confined" to a region of 
the phase space that is characterized by a NT. We still 
expect a convergence towards the equilibrium state pre- 
dicted in [5| eventually, but theprocess occurs on such 
long (unobservable) time scales 17j (further slowed down 
by coarsening processes) to make it unaccessible. It is 
nevertheless worth mentioning that if the IC contains one 
or more sufficiently-high breathers, the amount of energy 
that they carry is effectively frozen, while the structure of 
the metastable state depends on the remaining "excess" 
energy that is still free to diffuse between background and 
smaller breathers. Altogether, the scenario is quite rem- 
iniscent of what happens in spin glasses although a more 
quantitative analysis in necessary to clarify for example 
possible aging phenomena. 

Generation of NT states. NT states can be generated 
by moving in parameter space from an initial state in 
a region where equilibrium is characterized by a homo- 
geneous regime, i.e. positive temperatures. Two par- 
ticularly simple mechanisms to generate NT states are: 
localized dissipations with removal of mass (and energy) 
at the extrema of the chain [ll| or a free expansion in 
lattices of larger size. In the inset of Fig. [TJ we show an 
average trajectory that results from the first method of 
boundary dissipation in a chain with N = 200 sites. The 
free expansion method from a positive-temperature equi- 
librium state to a NT sate is demonstrated in the upper 
dot-dashed line in Fig. [5})). From Fig. [1] one can easily 
realize that a free expansion amounts to decreasing a un- 
til the T — oo line is crossed and the NT region is eventu- 
ally reached. Thus, in either BECs in optical lattices, or 
arrays of optical waveguides, this procedure amounts to 
prepare a standard equilibrium state at T > and leave 
it to spread in sufficiently larger lattice structures. This 
method is much simpler than the experimental schemes 
described in [12 . 

Conclusions. We have explored the dynamical behav- 



ior of the DNLS equation in the parameter region cor- 
responding to NT states. Numerical simulations show 
that for sufficiently uniform initial conditions, the sys- 
tem converges towards a stationary although metastable 
state characterized by a negative temperature and a fi- 
nite density of breathers. The asymptotic convergence 
to the single-breather equilibrium state occurs on much 
longer (and practically unattainable) time scales. We en- 
visage a quantitative theory that relates the time scales 
to the breather amplitudes, their stationary density and 
the system size so as to establish possible connections 
with, e.g. aging phenomena. 
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